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Abstract. Laser trapped nanoparticles have been recently used as model systems 
to study fundamental relations holding far from equilibrium. Here we study, both 
experimentally and theoretically, a nanoscale silica sphere levitated by a laser in a 
low density gas. The center of mass motion of the particle is subjected, at the same 
time, to feedback cooling and a parametric modulation driving the system into a non¬ 
equilibrium steady state. Based on the Langevin equation of motion of the particle, 
we derive an analytical expression for the energy distribution of this steady state 
showing that the average and variance of the energy distribution can be controlled 
separately by appropriate choice of the friction, cooling and modulation parameters. 
Energy distributions determined in computer simulations and measured in a laboratory 
experiment agree well with the analytical predictions. We analyse the particle motion 
also in terms of the quadratures and find thermal squeezing depending on the degree 
of detuning. 
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1. Introduction 

In a macroscopic system, thermodynamic quantities such as the work carried out during 
a thermodynamic transformation or the heat exchanged with a heat bath have well 
defined values due to the statistics of large numbers. For instance, if we repeatedly 
carry out a certain thermodynamic transformation always starting from the same initial 
state and following the same protocol, the work performed on the system will always 
be the same. In small systems, on the other hand, thermodynamic quantities typically 
fluctuate. Then the work and heat of a thermodynamic transformation, carried out, 
for instance, by stretching a single biomolecule in solution, need to be characterised 
with a statistical distribution rather than a single value. Even small systems, however, 
are subject to the basic laws of thermodynamics and, on the average, obey the second 
law usually formulated in terms of inequalities. As realised by Jarzynski, more specific 
results can be derived for the fluctuations of work and other quantities that transform 
the inequalities of thermodynamics into equalities [1, 2, 3, 4], which remain valid 
arbitrarily far from equilibrium. Such so-called fluctuation theorems have now been 
derived for several quantities, such as heat, work and entropy [5, 6], shedding new light 
on the significance of irreversibility and the second law at the nanoscale [7, 8]. Besides 
their fundamental importance, fluctuation theorems also provide the basis for the 
interpretation of single-molecule experiments [9, 10, 11] as well as for the development 
of novel non-equilibrium computer simulation methods [12], 

Experimentally, fluctuation relations have been studied in a variety of systems 
mainly in the over-damped regime, such as a particle dragged through a liquid [13] or a 
biomolecule in solution [10], where the system is strongly coupled to a thermalising 
environment. Recently, several experimental setups for the investigation of non¬ 
equilibrium fluctuations under low-damping conditions were proposed [14, 15, 16, 17]. 
Due to their weak coupling to the heat bath, such systems hold the promise to enable 
investigation of the statistics of non-equilibrium fluctuations in the quantum regime. 
Also, the precise control over the dynamics that can be achieved in such systems permits 
to construct situations in which microscopic reversibility does not hold. 

Here, we study, using theory, simulation and experiment a levitated nanoparticle in 
the low-friction regime [18]. In particular, we derive analytical expressions for the energy 
and phase-space distribution of the system in non-equilibrium steady states. Based 
on these distributions one can relate heat, entropy and energy to each other, thereby 
providing additional insight into the physics underlying the fluctuation theorems. The 
particle, consisting of a dielectric material, oscillates in a laser trap and is surrounded by 
a low-density gas, which exerts frictional and random thermal forces on the particle. The 
amount of friction can be controlled by changing the pressure of the gas. In addition, 
the particle is subjected to a nonlinear feedback cooling mechanism and a parametric 
modulation. Together, these effects allow to bring the oscillating particle into a variety of 
non-equilibrium steady states with tuneable parameters, turning such nano-mechanical 
oscillators into ideal test-systems for studies of stochastic thermodynamics. Based on a 
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Langevin equation written for the oscillating particle, we derive analytical expressions 
for the energy distribution in the stationary states and find that, under appropriate 
circumstances, our theoretical predictions agree very well with the energy distributions 
observed in the simulations. In addition, we find that in our experiments parameter 
fluctuations dominate the noise contribution from Brownian motion, which leads to and 
additional broadening of the experimental distributions. 

In addition to the levitated nanoparticle considered here, our model applies to 
other nonlinear oscillators, including ultra high-Q nano-mechanical oscillators fabricated 
from silicon nitride [19] and carbon nanotubes and graphene resonators [20]. The 
latter naturally exhibit nonlinear damping that is formally identical to our feedback 
mechanism. Thus, in addition to providing insights into thermodynamics on the 
nanoscale, the work presented here provides insight into the interaction of noise with 
inherent nonlinearities of nano-mechanical oscillators and the resulting amplitude and 
phase noise. Most notably phase noise, despite being an active topic of research for 
many decades, is still a pertinent topic today [21, 22, 23], since it plays a prominent role 
for the application of such systems as sensors and in timing and frequency control. 

The remainder of the article is organised as follows. In Sec. 2 we lay out the 
theory for the energy distribution of a nano-mechanical oscillator subject to friction, 
nonlinear feedback cooling and parametric modulation. Computer simulations are then 
used, in Sec. 3, to verify the theoretical predictions and probe the limits of the theory. 
In Sec. 4 we first describe our experimental setup and explain how we determine the 
relevant system parameters. We then present energy and phase distributions and discuss 
how they compare with theory and simulations. Some conclusions and an outlook are 
provided in Sec. 5. 

2. Theory 

2.1. Equation of motion 

We consider a particle of mass m oscillating in a trap with a Duffing potential 

V(q) = ^kq 2 + ^fkq 4 , (1) 

where q specifies the position of the particle, k is the trap stiffness, and £ is the Duffing 
parameter, which quantifies how strongly the trap deviates from a purely harmonic 
potential. Using the frequency D 0 = \Jk/m of the harmonic case, the total energy of 
the oscillator is given by 

E(q,p) = ^ mn 2 0 q 2 + ^ro^g 4 + (2) 

where p = mq is the momentum of the particle. The force due to the trap is hence given 
by 

^trap = -mQ 2 0 q - £mD[jg 3 . 


( 3 ) 
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Since the particle is immersed in a low density gas of temperature T, it experiences also 
a frictional force 


-^friction l^P (4) 

and the related fluctuating random force 

^random = \/2mT k B T w(t ), (5) 

where T is the friction constant, k B is the Boltzmann constant and w{t) is white noise. 
A feedback of strength rp acting on the particle with force 

-^feedback Pi (6) 

used to control the effective temperature of the center of mass motion of the particle 
and cool it far below the gas temperature T [18]. In addition, the particle is driven 
parametrically by periodically modulating the trap stiffness with frequency Vt m leading 
to the force 


-^drive f flli COS 


(7) 


where the modulation depth £ determines the intensity of the parametric driving. Taken 
together, these forces yield the following stochastic equations of motion for the motion 
of the particle in the trap, 

d q = —d t, (8) 

m 

dp = [— mVt^q — £mQlq 3 — Tp — Q 0 r]q 2 p + C m ^o cos(fi m t)g] di 

+ ^2mYk B TdW. (9) 

Here, W(t) is the Wiener process with 

(W(t)) = 0, (10) 

( W{t)W{t')) = min (t,t'). (11) 

Note that (W 2 (t)) = t for any time t > 0 and, thus, for an infinitesimal time interval 
d t one has ((dW) 2 ) = dt. The white noise w(t) appearing in the random force can be 
viewed as the time derivative of the Wiener process, w(t) = dW{t)/dt. 

In order to determine the energy distribution of the oscillator in the steady state, 
we now examine the time evolution of the energy generated by the stochastic equations 
of motion. To avoid multiplicative noise, i.e., a noise term with an amplitude depending 
on the current value of the energy, we consider the square root of the energy rather than 
the energy itself, 

e(q,p) = s/E(q,p). (12) 


Applying Ito’s formula [24] for the change of variables to e(q,p) we find that the change 
de during a short time interval is given by 


de 


’ pF(q,P,t) Tk B T f 

2 me 2e \ 



dt + \J 2 mT k B T-^—dW, 

2 me 


(13) 
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F(q, p , t) = —Tp — Q 0 r]q 2 p + (mQl cos(D m t)g 


(14) 


is the sum of the non-conservative forces consisting of the frictional force Faction, the 
feedback force Feedback and the driving force Fd r i ve - Note that the conservative forces, 
including the force due to the non-linear Duffing term in the energy, do not contribute 
to the energy change. 

The stochastic equation of motion for e, Equ. (13), explicitly depends on the 
position and momentum of the particle. To eliminate this dependence and obtain a 
closed equation depending only on e, we observe that the particle settles into a periodic 
motion with a frequency D that is not necessarily equal to the frequency Do of the 
unperturbed oscillator. Integrating Equ. (13) over one oscillation period r = 27r/D we 
we obtain the change Ae = fj de of e during the time r, 


Ae = - T 


p 


2 me 


dt — D 0 ?7 


q 2 p 2 

2 me 


dt + T/cb T 


2e 1 _ 


P 


2 me 2 ") 


dt 


+ f cos ^P dt + [ T / 


2 me 


2 me 


-dW. 


(15) 


To compute the integrals, we assume that during this time, which at low friction is 
short compared to the time for energy relaxation, the particle performs an undisturbed 
harmonic oscillation evolving according to 

q(t) = R cos (Dt + 0) p{t) = —mDF sin(Dt + 0), (16) 

where the amplitude R of the oscillation is related to e by R — \J2 /m(e/D). The 
phase 0 accounts for a possible phase shift with respect to the driving force, which is 
proportional to cos(D m t). Note that the oscillation frequency D is not necessarily the 
same as the frequency Do of the unperturbed harmonic oscillator. 

The central assumption, which allows to treat the motion of the system as that of 
an undisturbed oscillator during one oscillation period and eliminate the dependence 
on the rate of energy change on the phase space variables q and p by integration, is 
that the system evolves at nearly constant energy during one oscillation period. This 
condition is met if there is a separation of time scales between the time scale of the 
oscillation and the time scale for energy loss/gain. In other words, the relative change 
in energy A E/E occurring during one oscillation period should be much smaller than 
unity. The stochastic differential equation derived below for the time evolution of the 
energy, Equ. (27), provides a way to estimate for which ranges of the parameters T, 77 
and £ this condition holds. Analyzing each term on the right hand side of Equ. (27) 
individually, we find that the separation of time scale requires that T/D < 1, ( < 1 
and qkBT e s /mD 2 <C 1, where T e ® is the effective temperature of the oscillator. 

Carrying out the integrals over t, the first three terms in Equ. (15) yield 


a / re 

Ae =-r — 

2 4mD 2 


r]e 3 n 0 Tk B T 
r -|- ; —r 


4e 


( 17 ) 
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The change in e resulting from the driving (fourth term in Equ. (15)) is given by 

eCflo sin ( 7r % L ) [cos(20) sin(7rl|p) - (^) sin(20) cos(7r% L )] 


Ae" = 


SM 4-§- 


t.(1S) 


This expression is independent of time only if after one oscillation period the relative 
phase of the oscillation with respect to the periodic driving force is the same as at the 
beginning of the period. For the parameters studied here and a modulation frequency 
of 12 m ~ 212o, the oscillator locks to the modulation and oscillates with 12 = 12 m /2. We 
limit our considerations to this case in the following. Carrying out the limit 12 m —> 20 
in the above equation, one finds 

6 <(T2 q sin(20) 


Ae" = 


412 


-r. 


Finally, the last term in Equ. (15), 
Ae'" = ^2mTk B T 0 [ 


P 


2 me 


d W, 


(19) 


( 20 ) 


is a stochastic integral due to the noise term in the equations of motion. As a weighted 
sum of Gaussian random variables, Ae'" is also a Gaussian random variable with mean 

P 


and variance 


(Ae"') = y / 2mTk B T 


{(Ae ) z ) = 2mTk B T 


= 2mTk B T 


2 me 


(dW) = 0 


( 21 ) 



P(t)p(t) 


(dWdW) 


o Jo 


4 m 2 e 2 
p 2 r k B r 

dt = —-—r. 


4m 2 e 2 


Thus, the random variable Ae"' can be written in terms of the Wiener process as 
Ae'" = 


Tk B T 


W(t) 


( 22 ) 


(23) 


Putting things together, one obtains 

Te 


Ae = 


r?e 3 12 o T k B T 


e(Ql sin (20) 


r + 


Tk B T 


W(t). (24) 


2 4 mti 2 ' 4e 412 

Since the oscillation period r is short compared to all the time scale on which the energy 
changes, one can finally write the following stochastic differential equation for the square 
root of the energy £ 

Te 77l2 0 e 3 

’T - 


de = 


+ 


Tk B T eC!2o sin (20) 


dt + 


r k B T 


dW 


(25) 


4ml2 2 4e 412 

The corresponding Fokker-Planck equation [25] governing the time evolution of the 
probability density function P e (e,t ) is given by 
d P £ (e,t) d 


dt 


de 


Te 77 l2 0 e 3 
2 4ml2 2 


Tk B T eC12oSin(20) 


4e 


412 


PeY,t) 


Tk B T d 2 
4 de 2 


p,(c,t). 


(26) 
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In writing these two equation we have implicitly assumed that the phase <f between the 
modulation and the particle oscillation is fixed (or at least that it changes only very 
slowly in time). As we will show below, this condition is met very well particularly 
at low friction. Equation (25) implies that the time evolution of £ can be viewed as a 
Brownian motion in the high friction limit under the influence of an external force. Note 
that due to the integration over one oscillation period, this equation has e as its only 
time dependent variable while the dependence on other variables has been removed. In 
the following section we will use this equation to determine the energy distribution as 
well as the phase space distribution of the steady state generated by the parametric 
modulation and the feedback mechanism. 

Changing variables from e to E = e 2 and applying Ito’s formula [24] yields the 
corresponding stochastic differential equation for the energy, 

rjO 0 E 2 sin (20) 


d E = 


-r (e - k B T ) - 


dt+ ^/2EYk B TdW.{27) 


2 mO 2 20 

In contrast do stochastic equation of motion for e, here the noise is multiplicative, i.e., 
its amplitude is energy dependent. The corresponding Fokker-Planck equation for the 
probability density function P E (E,t) is given by 
dP E (E,t) d 


dt 


dE 


r(E-k B T) + ^ + Eca ° sia{2,p) 


2 mO 2 


20 


Pe(E,1) 


d 2 


TkBTg^EPEE.t). 


(28) 


2.2. Energy distribution 


The stochastic differential equation (25) has the form of the equation of motion 
describing the time evolution of a one-dimensional Brownian particle under the external 
force f(x) with large friction v at temperature T, 


dx = —f(x)dt + 


2 k B T 


-dW, 


(29) 


v V v 

where x is the position of the Brownian particle. The motion resulting from this equation 
of motion is known to sample the Boltzmann-Gibbs distribution 


Px(x) oc exp { -f3U(x )} , (30) 

where /3 = l/k B T is the reciprocal temperature and U(x) is the potential corresponding 
to the external force, f(x) = —dU/dx. 

By virtue of this isomorphism with over-damped Brownian motion, established by 
setting v = 4/T and identifying e with x, the determination of the energy in the non¬ 
equilibrium steady state of the driven oscillator turns into an equilibrium problem. One 
can then immediately infer that Equ. (25) samples the distribution 

P e (e) oc exp {—(3U (e)} , (31) 


where the potential 

U(e) = e 2 + 


r]O 0 e 4 
4 mTO 2 


— k B T In e + 


e 2 COo sin(20) 


2TO 


(32) 
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generates the force 


f( c ) = d 4 e ) = _ 2e - ?7 ^° e3 + kjiT sin(20) 


de ” mm 2 ' e m 

acting on the variable e. As a result, the systems samples the e-distribution 

>2 r/Oo 


P e (e) oc eexp < —/3 


J + W) Sio(20) ■ + 


( 33 ) 


(34) 


2rn /' 1 4mrsi 2 

Note that a small friction T corresponds to large friction v determining the time evolution 
of e and, thus, the energy E of the oscillator. By a change of variables from e to E, we 
finally obtain the probability density function of the energy E, 


Pe{E) = - exp 


-P 


Cflpsin(20) \ r]Q 0 

2m J 4 mm 2 


E 2 


The normalisation factor Z — J Pe{E)&E is given by 

limiTO, 2 , / I/3mm 2 (^ i £ngsin(2(/>) 


Z = 


Pffi 0 


-h 


rfil Q 


1 + 


2 m 


(35) 


(36) 


where the function h(x) is defined as 

h{x) = exp(x 2 )erfc(a:) (37) 

and erfc(x) is the complementary error function. Thus, the energy distribution is that 
of an equilibrium system with effective energy 


H = 


1 + 


COq sin (20) 


2TQ 


E + 


r]Q 0 


:E 2 


(38) 


4mm 2 

and configurational partition function Z. While the term proportional to E 2 is caused 
by the feedback cooling, the term proportional to E is affected only by the parametric 
modulation. 

According to Equ. (35), the energy distribution is Gaussian with a cutoff at E — 0. 
The maximum of the Gaussian is located at 


E = -- 


2 mm 2 
rf-lo 


1 + 


(tf 2 sin(20) 

2 m 


while its variance (neglecting the cutoff) is given by 


4 = 


2mm 2 k B T 

r]Q 0 


(39) 


(40) 


Hence, the width of the Gaussian does neither depend on the driving parameters nor 
on the phase 0. 
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Since for low friction the energy of the oscillator changes slowly, one can also obtain 
the full phase space density P qp (q,p ) from the energy density Pe(E). To determine the 
phase space density P gp (q,p), we consider the micro-canonical phase space distribution 
Pmc(q,P ; E) of the oscillator evolving at a given constant total energy E , 


Pmc{qi Pi E) 


1 

W) 


5(E(q,p)-E], 


(41) 


where 5(x) is the Dirac delta function and we have denoted the fixed value of the energy 
with E to distinguish it from the energy function E(q,p), which depends on the position 
q and the momentum p. The normalising factor g(E) is the micro-canonical density of 
states, 


9(E) 


= / dqdp 5[E(q,p) — E], 


(42) 


The phase space distribution of Equ. (41) would be observed for an oscillator evolving 
freely in the absence of feedback and without coupling to a heat bath. Since for the 
parameter ranges studied here the energy is essentially constant over many oscillation 
periods, the total phase space density P qp (q,p ) can be written by averaging the 
microcanonical distribution over the energy distribution, 


Pqp(q,p)= / d EP E (E)P mc (q : p-E) = 


d EtNP-S{E(q,p) - E}. (43) 

9(E) 


This linear superposition of micro canonical distributions is valid as long as the energy 
changes slowly on the time scale of the oscillation period. For the low friction constants 
and the small feedback strength studied here this assumption is met even under non¬ 
equilibrium conditions. Carrying out the integral yields 

P ^ ] = TOT < 44 > 

As further approximation, we now use the density of states g(E) = 27t/D 0 for the 
harmonic oscillator, thus neglecting the Duffing term of the potential in this part of the 
calculation, and obtain 


P«p(q,p) = ^Pb[s(?,p)]. 


(45) 


Inserting the energy distribution from Equ. (35) into this equation we finally find the 
phase distribution function 


P qp (q,p) = 


Do 


exp 


-P 


1 + cnfsin m\ E(q , p) + T n^ E{q , p? 


2TD 


4 mTD 2 


Note, however, that while we have neglected the Duffing term in the expression for 
the density of states, it is included in the energy appearing in the argument of the 
exponential on the right hand side of the above equation. 
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From the phase space density P qp (q } p ) one can obtain the distribution P q (q) of the 
position by integration over the momenta, 



ApP„(q, p). 


(47) 


In the absence of parametric modulation (£ = 0), one finds by carrying out the integral 

2 n 


a 


o / q 


p M) « \ 2 + 7]— — + £— exp 




/3mfl 0 

8 Tj 


a 


0 , -o q . ,q 

2 + r/— ( — + £— 


x K i 

4 


( fdmVtQ 

8p 


u 


o q 


2 + ^ 7 + ^7 


(48) 


where Ad /4 is a generalised Bessel function of the second kind. For simplicity, we have 
considered the case Oo = 0 here. A similar expression can also be derived for the 
momentum distribution. 


2-4- Relative entropy change 


As shown recently, a fluctuation theorem holds for the relative entropy change AS for 
a system relaxing towards equilibrium starting from the non-equilibrium steady state 
prepared by feedback cooling and parametric driving[15]. In this process, the feedback 
and the driving are turned off during the relaxation such that the system evolves freely 
and the dynamics is microscopically reversible. The relative entropy change AS is 
defined as the logarithmic ratio of the probability P[u(t)] to observe a certain trajectory 
u(t) and the probability P[u*(t)} of the time reversed trajectory u*(t), 

P[u(t)] 


AS = I 11 


(49) 


P[u*(t)]' 

Here, u(t) denotes an entire trajectory of length t including position and momentum of 
the oscillator and u*(t) denotes the trajectory that consist of the same states visited in 
reverse order with inverted momenta. Since during the relaxation detailed balance is 
obeyed, for the quantity AS a detailed fluctuation can be proven, 


P t (-AS)/P t (AS) = exp(—A S), (50) 

where P t (AS) is the probability density to observe the value AS at time t as 
determined over many repetitions of the relaxation experiment. For the relaxation 
process considered in Ref. [15] the relative entropy change is given by 

AS = PQ h + A0, (51) 

where Qh = — [Et ~ Eo] is the energy absorbed by the bath during the relaxation, and 
Eq and E t are the energy of the oscillator at time 0 and t, respectively. The quantity 
4>(q,p) is defined as as the logarithm of the stationary phase space distribution 

<i>(Q,p) = -fo.Pqp{Q,p) (52) 

and A(f is the difference of <f> at the beginning and the end of the trajectory, 


Ao - o, - (po¬ 


rn 
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Hence, the relative entropy change AS depends on the state of the system at the 
beginning and the end of the trajectory. 

In general, the steady distribution P gp (q,p) necessary to compute A0 is unknown. 
However, from the distribution derived for our model, Equ. (46), we find that for the 
relaxation from a non-equilibrium steady state generated by nonlinear feedback and 
parametric modulation, the relative entropy change is given by 

CHq sin(20) , a ?A) 


AS = -p- 


2m 


[E t ~ E 0 \ — P 


4 mm 2 


[E 2 t ~ Eg] , 


(54) 


Thus, our stochastic model allows us to express the relative entropy change during 
a relaxation trajectory in terms of the energy at the beginning and the end of that 
trajectory. Note that since no work is performed on the system, the heat Qh exchanged 
along a trajectory equals the energy lost by the system. Thus, in the absence of nonlinear 
feedback cooling, the relative entropy change is proportional to the heat and resembles 
relaxation form a thermal bath with effective temperature T cff = T/(l - c y ). By 
choosing parameters, one can therefore switch from a purely thermal situation with the 
phase space distribution of a harmonic oscillator (but with changed temperature) to 
a truly non-equilibrium steady-state with non-linear effects controlled by the feedback 
parameter q. 


2.5. Quadratures 

Parametrically driven nano-mechanical oscillators have been shown to support classical 
squeezed states in which the amplitude of the vibration in one phase is reduced with 
respect to the thermal equilibrium amplitude. To probe our oscillator for squeezed states 
we analyse its motion in terms of the so-called quadratures. For the oscillator driven 
by the parametric modulation Thrive = C m ^o cos(fi m t)g, we write the time evolution of 
the oscillator position as 

q(t) = R(t ) cos [fit + </>(t)], (55) 

where fi is the frequency of the particle oscillating at half the frequency of the driving, 
fi = fi m /2. Here, R(t) and <f>(t) are the amplitude and the phase of the particle, 
respectively, and the phase is measured with respect to the driving signal. Using the 
addition theorem for the sine-function, Equ. (55) can be written as the sum of two 
contributions, one in-phase with the driving signal and one out-of-phase, 

q(t) = R(t) cos <f>(t) cos (fit) — R(t) sin <f>(t) sin(fit) 

= X(t) cos(fit) — Y (t) sin (fit), (56) 

where the second line defines the in-phase component X(t) = R(t) cos <f>{t) and the 
quadrature Y(t ) = R(t) sin Together, X and Y are referred to as the quadratures. 
The quadratures can be computed from the time evolution of the position q(t) and the 
momentum p(t). The momentum of the particle is given by: 

p(t)/mQ— — R(t) sin[fit + (j)(t)] 

= — X(t) sin (fit) — Y (t) cos(fit), 


(57) 
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where we neglected the time derivatives of the amplitude and phase, since for an 
oscillator at low friction both the amplitude and the phase vary slowly in time. 
Combining this equation with Equ. (56) yields 

X(t) = q(t) cos (fit) — sin(Ql), 

mil 


Y(t ) = — qit) sin(f2i) — cos(f2l). (58) 

mil 

corresponding to transformation to a coordinate system that rotates clockwise with 
frequency with respect to the (q,p/mQ)- plane [26, 27]. In this coordinate system, a 
sinusoidal oscillation of frequency 11 is represented by a static point. 

Note that the amplitude and phase can be expressed in terms of the quadratures, 


R= VX 2 + Y 2 , 

(j) = arctan(iyx), (59) 


and that 

mil 2 /? 2 mil 2 
2 “ 2 

is the energy of a harmonic oscillator with frequency 12. 


(X 2 + Y 2 ) 


(60) 


3. Simulations 

In this section we verify the analytical expressions for the distributions of energy and 
positions by comparing them with simulation results. The simulations were performed 
for parameter values close to those of the experiments, which we will present and discuss 
subsequently. 

3.1. Simulation methods 

In our simulations, we integrated the Langevin equation of motion with the OVRVO 
algorithm of Sivak, Chodera and Crooks [28], which can be viewed as a stochastic 
generalisation of the velocity Verlet algorithm for deterministic dynamics [29]. This 
discrete time integration scheme uses a time step rescaling in the deterministic update 
step for positions and momenta to satisfy a number of desiderata proposed in the 
literature for stochastic integrators [30]. In all simulations we used a time step of 
At = 0.01 in reduced units. This time step is about 1/628 of the oscillation period. 
Test runs carried out with smaller time steps (At = 0.001) yielded identical results up to 
statistical errors. In most cases, the total simulation time was t = 10 7 corresponding to 
about 3 x 10 6 modulation cycles. For some parameters we carried out longer simulations 
of up to 3 x 10 10 steps corresponding to a total simulation time of t = 3 x 10 8 . All 
simulations were carried out for k^T = 1 , m = 1 , and k — 1. 

To facilitate comparison of the results of theory/simulation and experiments, in 
the following we use the thermal energy £ = ksT, the inverse frequency T = 1/flo and 
the particle mass M. = m as our basic units of energy, time and mass, respectively. 
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Accordingly, distances are measured in units of C = (1/flo) \fk^Tjm and velocities in 
units of V = \Jk^T/m. Hence, the unit of length is given by the variance of the position 
of the harmonic oscillator, ( q 2 ) = k^T/mOf = C 2 and the unit of energy is the average 
energy of the harmonic oscillator ( E) = k^T = £. The friction constant is given in 
units of such that it equals the inverse of the quality factor, Q = VL q /Y = 1/TT. The 
feedback strength r) and the Duffing coefficient £ have the dimension of 1/area and are 
measured in units of 1 /C 2 . The modulation depth £ is dimensionless. In the following, 
we use reduced units in which £ = T — M. — 1. 

3.2. Oscillator with feedback cooling but without parametric modulation 




Figure 1. Left: Energy distributions for different feedback strengths g without 
parametric modulation (£ = 0) for E = 0.0001, and £ = —0.022. The symbols are 
simulation results and the lines predictions according to Equ. (61). Right: Position 
distributions for the same parameters. The symbols are simulation results and the 
lines are theoretical predictions according to Equ. (48). 

We first consider the oscillator without parametric modulation (£ = 0.0) but 
subjected to feedback cooling. Without driving, the phase <f is not a relevant parameter 
and the expression for the energy distribution simplifies considerably, 

P £ (B)ocexp{-^B + i ^B 2 )}. (61) 

where we have assumed that the particle oscillates with fl = fl 0 . The first term in 
the exponential is the same as that of the uncooled oscillator, but the second term 
proportional to E 2 is due to the feedback loop and strongly penalises high energy states 
thereby cooling the system. The cooling effect is stronger for weak friction T and small 
frequencies Do- Several energy distributions obtained from simulations together with 
the corresponding predictions of Equ. (61) are shown in the left panel of Fig. 1. The 
simulations were carried out for a friction of T = 0.0001 and a Duffing parameter of 
£ = —0.022. Without feedback, rj = 0, the energy distribution is exponential, but for 
rj > 0 the E 2 term caused by the feedback suppresses high energies leading to a parabolic 
shape of the distribution in the logarithmic representation. In all cases, the theoretical 
predictions agree very well with the simulation results. Positions distributions for the 
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same set of parameters are shown in the right panel of Fig. 1. While without feedback 
the position distribution is Gaussian, the feedback quenches large deviations leading 
to a narrowing of the distributions. Also in the case of the position distributions the 
agreement between theory and simulation is excellent. 

3.3. Oscillator with parametric modulation but without feedback cooling 




Figure 2. Left: Energy distributions for different modulation depths ( without 
feedback cooling (rj = 0) for T = 0.01, £ = 0, k^T = 1, m = 1, k = 1, and 
fl m = 2 flo- The symbols are simulation results and the lines predictions of the theory. 
The theoretical predictions have been scaled by a factor such that they agree with the 
numerical results at high energies. Right: Phase distributions for different modulation 
depths £ obtained from the same simulations. 

We next turn to the oscillator with parametric driving but without feedback cooling. 
In this case, the energy deposited in the system by the modulation is removed only by 
the coupling to the gas as quantified by the friction constant T. If the particle oscillation 
is locked to the driving with a fixed phase q !>, the resulting energy distribution following 
from Equ. (35) is expected to be exponential, 

P e (E) <x exp {-/? (l + g} , (62) 

where we have assumed that the modulation frequency is Vt m = 2O 0 . For a vanishing 
Duffing parameter £ = 0.0, i.e., for a perfectly harmonic trap, the phase is expected to 
be 0 = —7 t/ 4 in the absence of thermal fluctuations [31]. If this is the case, the decay 
constant of the exponential is /?(1 — £T2o/2r). Hence, the decay constant is positive 
only for £ < 2r/D 0 . If the modulation depth £ exceeds this limit, the friction cannot 
remove the energy pumped into the oscillator by the modulation such that the oscillator 
energy keeps growing preventing the system from settling in a steady state. We indeed 
fold in our simulations that for £ > 2r/f2 0 the energy continuously increases. For 
weak driving, on the other hand, the energy distribution is expected to be exponential 
with the decay constant predicted by Equ. (35). Several energy distributions for this 
case are shown in Fig. 2. Note that we performed these calculations for a relatively 
large friction constant of Y — 0.01, because for lower friction it takes exceedingly 
long to sample all relevant energies. For weak driving, £ = 0.001 (red symbols), the 
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energy distribution is exponential as predicted by the theory. The negative slope of this 
distribution in the logarithmic representation is, however, slightly too large. The reason 
for this discrepancy is that the oscillation does not lock to the parametric driving as can 
bee seen in the distribution of the phase shown in the right panel of Fig. 2. The theory 
developed above, on the other hand, assumes a fixed phase of <f> = —7t/ 4 (for £ = 0). For 
£ = 0.001, the phase distribution is essentially flat implying that there is no preferred 
phase. As a consequence, essentially no heating occurs and the energy distribution is 
indistinguishable from the equilibrium distribution (black symbols). As the strength of 
the parametric driving is increased, a pronounced phase relation between driving and 
oscillation develops and two distinct peaks appear in the phase distribution at equivalent 
positions, one at = —7 t/ 4 and one at <f = — 7t/4 + n. Since the phase relation is more 
pronounced at high energies, in this regime the energy distributions shown in the left 
panel of Fig. 2 converge to the form predicted by theory. In the figure, the theoretical 
distributions are indicated by lines with logarithmic slope of — fd(l — £T2 0 /2r). For low 
energies, the phase relation is lost and the energy distributions have the logarithmic 
slope of the equilibrium distribution. Thus, the energy injected into the system by the 
parametric driving results in a longer tail in the energy distribution where it has the 
right phase relationship with the oscillation. In contrast at low energies, the form of the 
distribution is essentially unchanged with respect to the equilibrium distribution. 

3-4- Oscillator with feedback cooling and parametric driving 




Figure 3. Distributions of the phase f for friction constants T = 0.00001, 0.0001 
and 0.01 and for different Duffing parameter £. The simulations were carried out for 
g = 0.022, C = 0.03 and Q rn = 2 fl 0 . 

Next, we consider the oscillator with parametric driving and feedback cooling. 
To understand the energy distributions for this case, we first take a closer look at 
the statistics of the phase £>. In the derivation of the analytical energy distribution, 
Equ. (35), we have assumed a fixed phase between the modulation and the particle 
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Figure 4. Left: Most probable phase £>max as a function of the Duffing parameter £ 
for the friction constants F = 0.00001 (black), T = 0.0001 (red) and T = 0.001 (blue). 
The simulations were carried out for g = 0.022, £ = 0.03 and Ll m = 2 Qq. The symbols 
are simulation results and the lines are results of secular perturbation theory. Right: 
Number of full turns the oscillation fell behind the driving during the total simulation 
time of t = 10 7 as a function of the Duffing parameter £ for the friction constants 
T = 0.00001 (black), T = 0.0001 (red) and T = 0.001 (blue). 


oscillation. In practice, however, the phase </> follows a statistical distribution with 
a position and width that depend on the parameters, particularly on the Duffing 
parameter £ and the friction constant T. Several distributions of the phase obtained 
from our simulations for T and £ are shown in Fig. 3. These simulations were carried 
out for a modulation depth of £ = 0.03 and and a feedback strength of rj — 0.022, 
because these values can be realised in experiments. For all parameters considered 
here, the phase distributions are strongly peaked at a particular phase. The peaks 
are narrow for small friction and small Duffing parameters and broaden for increasing 
friction and non-linearity. Note that the Duffing parameters considered here are negative 
because the non-linearity is due to the shape of the focal intensity distribution, which 
is approximately Gaussian [32], Without non-linearity, £ = 0.0, the peak is located 
at <f> — — 7t/4 for all values of the friction constant. As one turns on the non-linearity 
by making the Duffing parameter more negative, the peaks become broader and shift 
towards more negative values. 

A closer analysis of how the phase depends on the Duffing parameter is shown in 
Fig. 4. The left panel of the figure shows the positions of the maximum of the phase 
distribution, i.e., the most likely phase 0 max , as a function of the Duffing parameter £ for 
different friction constants T. As can be inferred from the figure, the most likely phase 
( Lm determined from the simulations (symbols) follows exactly the form predicted 
by secular perturbation theory [31] (solid lines). While this theory neglects thermal 
fluctuations and cannot predict the entire phase distribution, it yields an accurate 
location of the maximum. 

Due to the thermal fluctuations, which lead to a broadening of the phase 
distribution, the oscillator might entirely loose the lock with the driving modulation 
and regain it only after falling behind by one entire turn of 2n. For the lowest friction 
studied here this never happens during a simulation of total time t = 10', but for higher 
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Figure 5. Left: Energy distributions for different friction constants T, for £ = —0.022, 
y = 0.022, £ = 0.03 and f l m = 2 LIq- The symbols are simulation results and the lines 
predictions of the theory. Right: Energy distributions for different Duffing parameters 
£ for t F = 0.00001, y = 0.022, £ = 0.03 and f2 m = 2 Do- The symbols are simulation 
results and the lines predictions of the theory. 


frictions, and in particular for large Duffing parameters, the oscillation may fall behind 
the parametric modulation several times. The number of times this occurs in the course 
of the simulations is shown in the right panel of Fig. 4 for different friction constants 
as a function of £. 

We now compare the energy distribution determined in our simulations for the 
oscillator with parametric driving and feedback cooling with the theoretical prediction 
of Equ. (35). To do that, we identify the phase 0 occurring in the theoretical expression 
with the most likely phase </> max determined in the simulations. Energy distributions 
obtained for friction constants ranging from T = 10” 5 to T = 10~ 3 are shown in 
the left panel of Fig. 5. In all cases, the system was driven at Vt m = 20o and the 
Duffing parameter, the feedback strength and the modulation depth were £ = —0.022, 
r/ = 0.022, £ = 0.03, respectively. While for high friction the theoretical predictions 
deviate considerably from the energy distributions determined in the simulations, most 
likely due to the lack of a stable phase relation, very good agreement is obtained for low 
friction, where phase distributions are strongly peaked. This excellent correspondence 
is confirmed by the energy distributions shown along with theoretical predictions in the 
right panel of Fig. 5 for different Duffing parameters at low friction. Thus, the position 
and the width of the energy distribution in the non-equilibrium steady state generated 
by driving and cooling at the same time can indeed be controlled independently by an 
appropriate choice of parameters. 

3.5. Quadratures 

Finally, we take a look at the distribution of the quadratures X and Y for different 
driving frequencies. Scatter plots of the quadratures obtained at different driving 
frequencies and for different values of the friction constant are shown in Fig. 6. From 
left to right, the driving frequency kl m is slightly below 2D 0 ; equal to 2D 0 an d slightly 
above 2 Dq- As in previous simulations, the parameters were £ = —0.022, p = 0.022, 
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Figure 6. Scatter plot of the quadratures X and Y for the friction constants 
T = 0.00001 (black), 0.0001 (red) and 0.01 (blue) for = 1.98 (left), = 2fi 0 
(center) and f2 m = 2.02 (right)- The simulations were carried out for g = 0.022, 
£ = -0.022, and C = 0.03. 

and ( = 0.03. At Q m = 2fi and low friction the quadratures of the driven system 
are Gaussian with equal width along the two quadrature axes. Thus, they resemble a 
thermal state, albeit, displaced from the origin. In contrast, for driving frequencies off 
2O 0 , the distributions are deformed, indicating the occurrence of classical squeezing. 

4. Experiments 

In this section, we discuss how to retrieve the energy and phase of a trapped nanoparticle 
from discrete measurements of the particle positions. From the retrieved energies and 
phases we reconstruct the energy and phase distributions and compare them to the 
theory and simulation results presented in the previous sections. This allows us to 
extract the experimental parameters, which are detailed in Table 1. While the maxima 
of the distributions are in good agreement with our theory and simulations, the width 
of the experimental distributions is significantly broader due to parameter fluctuations 
not taken into account in the theoretical considerations. 

4-1. Experimental configuration 

In our experiments we use a silica nanoparticle trapped at the focus of a single beam 
optical tweezers. The optical tweezer is formed by a 1064 nm laser beam (~ 35mW) 
focused by a NA = 0.9 objective, which is mounted inside a vacuum chamber. The 
particle motion is recorded with an additional colinear laser (780nm) and three balanced 
photodetectors. A home-built electronic circuit is used to generate the feedback signal 
(rj), while a frequency generator serves as the parametric modulation signal (£). The 
approximately Gaussian shape of the optical potential is responsible for the trap 
anharmonicity (£) [32]. The detectors and the size of the nanoparticle are calibrated from 
measurements of the power spectral density of the particle motion at 5.1mBar. At this 
pressure the Q-factor is high enough to resolve the three spatial modes, while broadening 
effects due to nonlinear mode coupling are negligible [32], For further details of the 
experimental configuration and calibration procedure see Refs. [33, 34], Subsequent 
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Parameter 

Value (phys. units) 

Value (dimension less) 

error (%) 

a 

82 ± 4 nm 

2.3 x C 

4 

m 

5.2 ±0.7 x 10" 18 kg 

1 x M 

13 

V 

3.9 ± 1.3 /im~ 2 

4.9 x 10“ 3 x Cr 2 

34 

£ 

—5.4 ±1.1 pm~ 2 

-6.9 x 10~ 3 x C~ 2 

20 

r 

2n x 8.1 ±0.2 x 10" 3 Hz 

6.25 x 10~ 8 x T" 1 

3 

Q 

1.54 ±0.03 x 10 7 

1.54 x 10 7 

3 

Ho 

2vr x 125.12 ±0.05 kHz 

1 x T" 1 

0.04 

c 

16.1 ± 1.3 x 10” 3 

16.1 x 10" 3 

36 


Table 1 . Overview of experimental parameters. The second column lists the 
parameter in SI units with their respective experimental uncertainties, while the third 
column shows the experimental parameters in dimensionless units. For the scaling to 
dimensionless units see section 3.1. The last column lists the relative uncertainty of 
the experimental parameters. 


measurements are carried out at 1.2 x 10 _5 mBar. 

While our theoretical model is one-dimensional, the particle in the experiment 
moves in three dimensions along three main axes. The three axes are determined by 
the symmetry of the laser focus. However, there is no direct coupling between the three 
spatial modes. In addition, feedback cooling reduces the amplitude such that also the 
nonlinear coupling becomes very weak. Therefore, our one-dimensional model is a very 
good approximation for the particle motion along one of the three main axes. 


f.2. Amplitude and phase estimation 


The particle oscillation frequencies along the three main axes are well separated and 
dont overlap. Therefore, we can apply the maximum likelihood estimation for a single 
tone signal, that is a signal containing only one frequency component. The maximum 
likelihood estimation of the oscillation amplitude and phase of a single tone signal q{t) 
is given by [35] 


-Rml — |Aj(H)| (63) 

0ml = arg [exp(-iQ 0 t 0 )A q (Q)] (64) 


where fi/ 2tt is the estimated frequency of the signal, t 0 is the time origin and 


A q (uj) 


1 

N 


N—l 

q n exp(— icon At). 

n =0 


(65) 


is the discrete Fourier transform of q evaluated at oj. Here, q n = q(t n ) is the measurement 
sample of the time trace at time t n = to ± nAt, N is the number of samples entering 
the estimation and At is the sampling interval. The estimation of the amplitude and 
the phase relies on precise estimation of the frequency fi. We estimate fi by maximising 
(65) with respect to cv, i.e. H(fi) = max(H(u;)). The width of the function A{oj) 1 and 
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thereby our ability to localise the maximum, depends on the length of the time trace 
q{t). Therefore, we use a long time trace measured over T meas = 0.1s and sampled at 
625 kilosamples/second to estimate hi. Subsequently, we use that value of and Equs. 
(63) and (64) to estimate the instantaneous amplitude and phase from short parts of 
that same time trace. The short parts of the time trace contain N = 160 samples, 
corresponding to an integration over 32 particle oscillations. This constitutes a good 
compromise between sufficient data points for an accurate estimation of R and </>, and 
fast time resolution to resolve the dynamics of the energy and phase fluctuations. Note 
that maximising (65) allows us to estimate the frequency with much better accuracy 
than 1 /T meas .. 

The absolute phase of a harmonic oscillator is a time delay with respect to some 
time reference. Without such a time reference the absolute phase is arbitrary and has 
no meaning. However, the relative phase between two oscillators is meaningful, because 
one oscillator serves as a time reference to determine the phase of the other oscillator 
with respect the first oscillator. Formally, this is expressed as 

= {(Pp-^[(P m +k2n}), ( 66 ) 

where A p and A m are the Fourier transforms of the two signals, respectively (c.f. (65)), 
and k is an integer which takes into account that the phase is only determined up to 
modulo 27 t. Note that the exponent Q p /Q m takes care that ( 66 ) does not depend on 
to- Without loss of generality, we set <f> m = 0, i.e. we choose our time origin such 
that it coincides with a maximum of the signal with frequency f l m . For the special 
case of a parametrically driven particle, which oscillates at half the frequency of the 
parametric modulation (fi m = 2Q p ), we get A0 = <f> p — kn. Therefore, the above 
method allows to estimate the relative phase between the particle oscillation and the 
parametric modulation up to a multiple of i r. 

f.3. Parameter estimation 

We measure the distribution of the energy and phase for modulation at VL m /2n = 
247,248,249, and 250 kHz. Each distribution is obtained from 100 time traces of 0.1s 
duration. Fig. 7 shows the maximum values of the energy and phase distributions shown 
in Fig. 8 and a fit to secular perturbation theory [34, 31]. While independent fits to the 
energy and phase, shown in blue and red, respectively, yield excellent agreement with 
the theoretical model, we cannot fit a set of parameters that would agree with both the 
energy and the phase. Note that the phase fit includes a constant phase offset f> 0 = 50° 
to account for the finite response time of the intensity modulator and delays in the 
electronics. Averaging the results from the independent fits to energy and phase yields 
£ = —5.4 ±1.1 /rm~ 2 , 77 = 3.9 ± 1.3 /im ~ 2 and £ = 16.1 ± 5.7 x 10 -3 . The theoretical 
curve for the parameters obtained by the energy and phase is shown in green together 
with numerical simulations using the parameters summarised in Table 1. 


A 0 = arg 


A ri 


4 s 

r±r 
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Figure 7. Left: Most likely energy. Right: Most likely phase. The black and 
green circles are the experimental data points and simulation results, respectively. 
The blue and red solid lines are the theoretical predictions for parameters obtained 
from independent fits to the energy and phase, respectively and the green solid line is 
the theoretical prediction for the averaged parameters. 



Figure 8. Left: Experimental distributions of energy Right: Experimental 
distributions of phase. The circles are experimental data points and the solid lines 
are Gaussian fits. The maxima correspond to the data points shown in Fig. 7. 


The main uncertainty in the determination of the experimental parameters arises 
from the estimation of the particle mass and the resulting uncertainty in the voltage 
calibration and from parameter fluctuations, which we discuss in the next section. As an 
independent measurement, we also measure the energy distribution without parametric 
modulation (( = 0). A fit of the energy distribution to Eq. (61) yields r) = 4.5±0.9 pm -2 , 
in good agreement with the previously determined value. 

f.f. Distributions 

Fig. 8 shows the experimental energy and phase distributions fitted with a Gaussian. 
As predicted by our theory and simulations, the distributions are Gaussian and their 
widths depend only weakly on the modulation frequency. Fig. 9 shows the widths of the 
distributions obtained from the Gaussian fits in Fig. 8 and from numerical simulations. 
For comparison, we also show the theoretical prediction according to Equ. (40). The 
broadening of the distributions has two contributions, thermal motion and parameter 
fluctuations. 
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Figure 9. Left: Widths of energy distributions. Right: Widths of phase distributions. 
The black circles are experimental data points obtained form the Gaussian fits in Fig. 
8. The green circles are simulation results and the green solid line is the theoretical 
prediction Eqn. (40). Note that the experimental values are significantly larger than 
the theoretical ones and are scaled by a factor of 0.1 to fit them into the same plotting 
range. 


Thermal motion of the resonator, caused by residual air molecules, enters directly 
as a random white noise, which we considered in our theoretical model. In addition, it 
enters indirectly through amplitude-phase conversion [36]. The latter contribution has 
not been considered in our theoretical model but is naturally present in the numerical 
simulations. Amplitude-phase conversion refers to the interdependence of energy and 
phase (c.f. (39)). Therefore, fluctuations in the phase cause fluctuations in the energy 
and vice versa. This leads to a broadening of the distributions near the instability 
boundaries, where the deviation of the numerical simulation from our model is largest. 
Within this range, on the other hand, this interplay manifests itself as sidebands in the 
power spectral density of the particle position [34], 

In addition to Brownian motion, parameter fluctuations broaden the experimental 
distributions [37]. The experimental parameters fluctuate due to laser intensity and 
polarisation fluctuations and also due to the nonlinear coupling with the other two 
degrees of freedom, which were not considered in our model [32, 34]. Noise in the 
feedback electronics and modulator gives rise to further broadening. In general, 
broadening due to fluctuating parameters dominates broadening due to Brownian 
motion. As a consequence, the measured width of the energy and phase distributions 
a E = 78 ±3 x 10” 3 k E T and = 1.7 ±0.1 x 10” 3 7 t, respectively, are approximately one 
order of magnitude larger than the theoretical values 5.1 x 10 3 k E T and 0.15 x 10 3 7 t, 
averaged over the range of detunings of the experimental data. To identify the noise 
sources responsible for the deviation from theory one can deliberately introduce noise 
and systematically study its effect on the measurement outcome. 

5. Conclusion 

We have developed a stochastic model for the dynamics of the energy of a nonlinear 
nanomechanical oscillator subject to parametric modulation and nonlinear damping. 
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Under these conditions the oscillator attains a non-equilibrium steady state. Our model 
allows us to predict the energy distribution of the steady state. The steady state 
distribution is intimately related to fluctuation theorems, which describe the statistical 
properties of the system for transitions between different states [15]. Consequently, our 
work opens the door to test these fluctuation theorems in different scenarios. 

We confirmed the validity of the model by extensive numerical simulations and 
found excellent agreement with our theory. In addition, we performed experiments 
with a levitated nanoparticle. While the measured mean energy and phase are in close 
agreement with the numerical simulations, their distributions are broadened due to 
parameter fluctuations that are not accounted for in the theory and are subject to further 
investigation. Besides quantifying additional noise sources experimentally, future work 
includes the development of a more generalised model including a stochastic model for 
the phase and incorporating other white and non-white noise sources, resulting from 
fluctuating parameters [23]. 
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